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We derive low-temperature properties of the large-C/ Hubbard model in two and three dimen- 
sions from exact series-expansion results for high temperatures. Convergence problems and limited 
available information prevent a direct or Pade-type extrapolation. We propose a new method of 
extrapolation, which is restricted to large U and low hole densities, for which the problem can be 
mapped on that of a system of weakly interacting holes. In the new formulation an extrapolation 
down to T — can be obtained, but it can be trusted for the presently available series data for 

.^4- | fit < 20 and for hole densities nu ^ 0.2 only. Implications for the magnetic phase diagram are 

Qv^ ■ discussed. 
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I. INTRODUCTION 



The single-band Hubbard model is presumably the simplest model for describing the behavior of correlated electrons 
in a solid. Examples of its applications are its initial use to describe magnetism in transition metals, and, most recently, 
theories of high-temperature superconductivity. Unfortunately, while it seems likely that for the latter phenomenon 
more complex models are needed, even this simple model is not nearly well understood. For one dimension some 
rigorous results are known, but in higher dimensions the main results have been obtained from Monte Carlo and 
£■ — . finite-lattice calculations only. 

We are interested in deriving .magnetic properties for the case of large U on a square or simple-cubic lattice. A 
well known theorem by NagaokaEl states that a Hubbard model on a bipartite lattice with one hole and at infinite U 
has a ferromagnetic ground state. Many authors have investigated whether this one point in the phase diagram 
is part of a whole region of ferromagnetic behavior. Various methods are being used for this purpose!], including 
exacti diagonalization of small systems^rMonte Carlo simulations, mean-field, and cluster expansion methods. Two 
of usS as well as various other authorstffl have used the last method to calculate high-temperature series expansions 
for the square and simple-cubic lattices. Expressions have been obtained for various thermodynamical quantities, 
such as the free energy, the magnetization, the magnetic susceptibility, and also for the pair-correlation functions 
between the z components of the spin at specific sites. These expressions show very well convergent behavior for high 
O ■ temperatures (kT/t > 2). The aim was to find indications for the onset to ferromagnetic behavior at low temperatures 
^. \ by extrapolating the results of the series expansions. Indeed, these indications can be found, as is shown in Refs. y 
and 0. However, predictions for the ground state, based on these results, are highly unreliable. Due to the fact that 
we only have five terms in the series expansions (0th, 2nd, 4th, 6th, and 8th order terms), we found it impossible 
to rely on standard extrapolation methods like Pade approximants. The obtained r^ults for different extrapolations 
vary too much to be able to derive any reliable extrapolated value. Henderson et al.Q try to find an indication for the 
expected divergence in the uniform susceptibility 



*- = /3 ^ lntre ^ 



(1) 



h=0 



by looking for zeros of x fm . The character of the series expansion, shown for infinite U as a function of j3t in subsequent 
approximations in Fig. nl is such that xTm (fo) is likely to diverge very quickly for f3t > 1, to plus and minus infinity 
alternatingly. This means that zeros are to be found for j3t > 1 in the fourth- and eighth-order approximations, but 
no zeros exist in the second and sixth orders. 

We feel that there is no reason to believe that the fourth- and eighth-order results should be more reliable than the 
others. 

Furthermore, we have also constructed the antiferromagnetic susceptibility % af by including a staggered-field term 
in the Hamiltonian, and we have calculated its divergence in the same way as described above. In Fig. we compare 
the Curie temperature Tq as a function of the particle density, for various values of t/U, to the Neel temperature In, 
for calculations up to eighth order (Tq and In are defined by Xfm ( n i ^ ^c) = and X&1 ( n i U, In) = 0, respectively). 
As the Neel temperature is higher than the Curie temperature for the parameters shown, one should conclude that 



the system goes into an antiferromagnetic state before the ferromagnetic transition is reached. However, regarding 
the character of the series expansion, as illustrated in Fig. [l| it is clear that the plots can not be trusted qualitatively, 
let alone quantitatively. 

In this paper we will consider a method that does not encounter these problems of extrapolation to low temperatures. 
In this method the density of holes is used as a small parameter. The high-temperature results are expressed in terms 
of an effective density of states for holes (as was done before by Brinkman and RiceO) , and extended to interactions 
between hole levels. With this density of states expressions for the free energy of the thermodynamic system can be 
obtained in the whole range of temperatures. In Sec. y| we will define a partition function for the holes and express it 



in terms of an effective chemical potential for the holes. In Sec. Ill we derive the density of states for non-interacting 
holes, and we determine its moments, for infinite U. We present an improvement on the non-interacting hole picture 
in Sec. LVL where we consider interacting holes by introducing a Fermi-liquid-like interaction in energy space. In 
Sec. [V] we show how to use the density of states to calculate zeros of the inverse susceptibility. Sec. VI deals with the 



non-interacting hole approximation applied for finite U. In Sec. VII we show our conclusions for the magnetic phase 



diagram, and we discuss the method in Sec. VII] 



II. HOLE FORMULATION 

We consider the Hubbard Hamiltonian 

H = Hkin + TYlocal - /» X! n ™ ~ h Z^2 an i° ( 2 ) 

with 

Hkin = -t Y^ C L C j<? ( 3 ) 

i 

where t is the hopping integral between nearest neighbors, U denotes the on-site interaction strength, fi is the chemical 

potential, and h is the strength of an external magnetic field. The operator c ia (c ia ) creates (annihilates) a particle 

with spin a at site i, and n icr = c ia c ia counts the number of particles with spin a at site i. 

To investigate the thermodynamic properties we want to calculate the grand canonical partition function 

Z gr = tr e~ m (5) 

For a system consisting of N sites we can rewrite this as 

IN 

^=E e ^ s ^ (6) 

with Zn 3 the canonical partition function for N s = iVj + N^ particles: 

N a 

Z Ns = g em^-N,) £ e -/3 £ < N - N ^> (7) 

JV T =0 3 

Here {ea ' s ' } is the set of eigenvalues of Hkin + Wi OC ai for ^\ up spins and N^ down spins on N sites (note that 
the Sj are functions of t and U only). Via the grand potential 

Q = ~lnZ gt (8) 

we can now derive the other thermodynamic quantities by means of the usual manipulations, e.g. the particle density 

_ (N s ) 1 d(3Q 



N N dj3fi 



(9) 



or the susceptibility as given by (P. 

In order to approach to lower temperatures in the limit of strong interactions and near half filling, we are going to 
express the partition function in terms of an effective chemical potential for holes. We associate the kinetic part of 
the Hamiltonian with the motion of the (dilute) holes, and its magnetic part with the background of spins. Thus, we 
have to divide out the spin degrees of freedom to obtain the canonical partition function for the holes: 



Here we define the number of holes 



Z h Nh = Z N _ Nh e^ N - N ^ (10) 



N h = N-N s (11) 



and we introduce a parameter ehf which can be viewed as the free energy per spin in the absence of holes (i.e. at half 
filling; naturally, Zq = 1): 

Chf ~ ~~NB lnZN ^ 

= — - ln(2 cosh^/i) (13) 

P 

for infinite U . The grand canonical partition function for the holes then is 

Z\ T = Z gI e p(e "-^ N (14) 

= Y,z^y {ih! ~ ti)Nh (15) 

N h 

suggesting the definition of an effective chemical potential for the holes (Cf. Eq. (g)): 

Mh = £hf - A* (16) 

With this definition we can rewrite (O) as 

lnZ gr = -/?/i h iV + ln4 1 r (17) 



Note that expression (T13J) for £hf is exactly true only in the case of infinite U, as the interaction then prevents 
particles from occupying the same site. Note also that we do not define the number of holes as the number of sites 
where no particles are present (a definition which seems obvious), because the interpretation of Eqs. ( |10|) and ( |15| ) 
would then become problematic for finite U. However, if U is very large, as we assume, a pair of electrons located 
on the same site causes a very high energy, and the contribution of the corresponding hole to the kinetic part of the 
Hamiltonian is some orders of magnitude smaller than the contribution of a 'real' (non-removable) hole. Therefore we 
will use ( |ll|) also in the case of large, finite U, and we will show that this leads to terms to be added to the expressions 
for infinite U of order jj or higher. 

III. CONSTRUCTION OF THE DENSITY OF STATES FOR INFINITE U 

Let us consider a system near half filling, with, for simplicity, infinitely strong coupling U (the case of finite U will 



be treated in Sec. VI). We assume that the system can be described in terms of the kinetic energy of non- interacting 



dilute holes and the magnetic energy of the background particles. We define 

||U Jdep(e,(3h)e-^ (18) 

where p(e, j3h) is the spectral distribution of the energy levels of one hole in an otherwise half-filled system. We take 
p to be normalized to one. 

Although we said before that we divide out the magnetic degrees of freedom in the spin background, there is still 
a dependence of p on the magnetic field h. It is not easy to see how the hole motion depends on the field exactly, 
but one can easily understand why this dependence exists: a magnetic field influences the distribution of the spin 
background, which in turn determines the behavior of the hole. The hole motion depends on the field only indirectly, 



and the mechanism that governs the hole dynamics can in fact be much better described in terms of the average 
magnetization of the spin background than in terms of the field. It is important to understand that, in this picture, 
one has to treat the spin background as if it were at half filling, with the dilute holes subjected to its magnetization. 
Therefore we change variables at this level from /3h to the magnetization per spin to. This change is easily performed 
by a Legendre transformation 



£hf(wi) = Chi(0h) + mh (19) 



yielding 



where p(e, to) is obtained from p(e, (3h) via 



zl 

N 



dep{e,m)e- ,3t£ (20) 



to = -— (21) 

which becomes to = tanh(/3/i) for infinite U. 

With this definition we can write down a first approximation for the grand canonical partition function. A one-hole 
level can be occupied, with a Boltzmann weight e~@ t6 (where we write te to make the dependence of the hole energy 
on t explicit), or it can be unoccupied, in which case there is an electron in the system with Boltzmann weight 
e -PA% ( w ith the magnetic energy included in p^). Thus, in the case of non-interacting holes we have (dropping the 
TO-dependence of p) 

In Z gv = N f dep(e) hx{e~ 0te + e~ p ^ ) (22) 



or equivalently, using fll1 

]nZ* =N f dep{e) In (l + e ~P^-^)\ (23) 

This equation becomes exact in a one-dimensional system, as in that case the holes can not disturb the magnetic 
background of the particles, thus being really non-interacting, and also in a ferromagnetic system (at to = ±1), for 
similar reasons. In other, higher-dimensional systems ( |23|) is only correct to first order in e^ h . We make an expansion 
of the right-hand side with respect to the small parameter e 13 ^ to obtain 

In Z% = N I dep(e) (e^ h e~ me - I e 2/5 ^ e - 2/3t£ + . .\ (24) 

Comparing this to a similar expansion of the logarithm of Eq. ([La) we see that this is consistent with the definition 
of the density of states in the first-order term. 

Now, as an illustration of the calculation, let us have a look at the form that Eq. (||) actually takes when evaluating 
it for this system by means of a cluster expansion method. In Ref. a general formula is presented for finite U, which 
for infinite U reduces to 

, rr oo 2n-l m /-\(n) /3(m^+lh) 

^ = mz gr , + EW 2 " E E % oYn (25) 

iV „=1 m=0 l=-m ^^-0) 

Here i? gr ,o denotes the partition function for a system consisting of only one site: 

Z gr , = 1 + 2e^ cosh(/?/i) (26) 

(n) 

The yt m \ are coefficients, which can e.g. be calculated by means of a cluster expansion method. By substituting /ih 
for p, using ( jlq ) and (|l3|), and expanding in the small parameter e^ 11 , we can obtain an expression for the grand 
potential for the holes again, now in the form of a series expansion: 

In Z h °° °° 

_JE = J2 e^^E^*) 2 "^^"^) ( 27 ) 

N h =l n=0 



where 



(_l)(p-i) 
fi(p,<U) = ^ (28) 



and 

In — 1 m 



_ 1 \ (-l)PQ w e /^ 
n(p,n,fc)-X; E {p + m-2n) { -2coMPh)r ^ 

m=0 (=— to 

for n^O. Finally, we obtain a relation between the coefficients fl(l,n, h) and the moments of p(e) (= p(e,[3h)) by 
expanding (fi"q) in powers of (3t: 

-w = EC 9 *: ) nL rr M^K* (so) 

Thus, we see from Eq. (|3^) and the first-order term in ( |27| ) that we have 

dep(e)e 2n = (2n)!f7(l,n,/ l ) (31) 

for the even moments of p, all odd moments being zero. Although we have restricted ourselves to the case of infinite U 



here, this expression can easily be extended for the case of finite U, as we will see in Sec. VI . U then enters the equation 
as a parameter at the right-hand side. 

For infinite U there is another, faster way to calculate these moments. They can then be expressed directly,4n the 
number of possible paths in state space for a system with one hole. This has been done by Brinkman and RiceB, who 
calculated the first 10 moments of the density of states for ferromagnetic, antiferromagnetic and paramagnetic spin 
backgrounds on a simple cubic lattice, and by Meshkov and BerkovQ, who calculated 16 moments for the same spin 
backgrounds on a square lattice. In Appendix A, we outline a method which enables us to enumerate the paths in an 
efficient way, and by which we have extended the results for the square lattice to 22 moments. These moments are 
presented in Table |i| for m = and m = 1 , corresponding to a paramagnetic and ferromagnetic system, respectively. 
We have also extended the results for the simple cubic lattice, to 16 moments; they are available on request. 

We now approximate p{e) by a polynomial which we fit with the moments. In this way we calculate an approxi- 
mation for the density of states, to different orders, in order to get an impression of the convergence of subsequent 
approximations. In Figs. 0a and pb we show the result for a paramagnetic and a ferromagnetic system. 

For m = 1 the exact density of states is known: 

1 
2~^ 2 "W V4 



^) = o^(1-(t) ) ( 32 ) 



with K the complete elliptic integral of the first kind. It has an integrable singularity at e = 0. This is difficult to 
approximate and causes some oscillations away from e = 0. Convergence towards the exact result is rather good. For 
in = convergence is very good, as from 14th order on the difference between subsequent approximations becomes 
very small. 

Meshkov and Berkov fit the density of states by postulating that the integral of p 2 be minimal ('smoothness' crite- 
rion), using a discretized p. They claim that this method gives faster convergence then a polynomial fit. Comparing 
their results for the ferromagnetic density of states with the exact result and the results presented here, however, one 
may question that claim. We feel that the polynomial fits, when using an equal number of moments, give similar or 
even better results, which are also easier to handle in further calculations. 

Before calculating various quantities which can tell us something about the low-temperature properties of the 
system, we will in the next section consider a method to improve the approximation of the density of states by 
including interactions between the holes. 

IV. INTERACTING HOLES 

The crucial question is to see for which domain of hole densities the assumption of independent holes is justified. 
This range can be determined from an estimate of the interactions between the holes. Very similar to the theory of 



the classical dilute gasH, the interaction can be deduced from the two-hole partition function as defined by (|10|) for 
Nh = 2. It is a matter of choice how to represent the hole interaction. One could think of a spatial representation, but 
one must realize that in this strongly quantal system the interaction is non-local, which complicates the transparancy 
of the representation substantially. Having the one-hole system represented by a density of states it is natural here 
to choose an interaction between energy levels. First we formulate the interaction in terms of discrete levels and then 
we take the continuum limit as in Eq. ( |22|) . The discrete version of this expression can be written in terms of levels 
Si, distributed according to the density p(si): 

Zg) = e -^ N Y, e^E^-'"')"* (33) 

where {rii} with m — 0, 1 is the occupation of the levels e*. We have given the expression a superindcx 1 to indicate 
that Z gr matches Z gT up to the one- hole terms. The next approximation can be of the form 



.u w jvv--^(E 4 ( fe «-^+E (u) A 



Z$=e-to* N Y, e V ^> — >) (34) 

where fy accounts for the interaction between the levels Si and Ej. The second term in the exponent is a sum over 
all pairs of levels (i,j). In the energy space a distance between levels does not seem to be a measure for the strength 
of the interaction as in real space, where interactions usually decay sufficiently fast with the distance, such that the 
sum over pairs does not increase with the square of the number of elements, but only linearly as is necessary for a 



thermodynamic system. In order to make the exponent in (34) of the correct thermodynamic behavior the interaction 
should therefore decrease with the size of the system as 



t 

AT 



fv = —&J ( 35 ) 



with <fiij of order unity. An additional advantage of (pq ) is the fact that interactions of this type can be handled 
rigorously in the thermodynamic limit by the mean-field theoryQ. Thus we can write 

In Zg> = -fa h N + J2 ln (! + e _/J * ) + § E &M^M^) (36) 

1 («'.i) 



where the £ are the shifted energy levels 



£i = te t - n h + — y~] 4>i n(£j) (37) 






and n(e) is the fermi occupation number 



1 



n & = Tl-Si ( 38 ) 



(2) 

Now the interaction <pij must be chosen such that Z gr produces the correct two-hole partition function. Expanding 
Eq. (JM) with respect to the number of holes 



N h =J2m (39) 



and using Eq. (14) we find 



Z\ = V^ e -M£i+t£j + -k<t>ij) (40) 

In our high-temperature expansion we have no direct information on Z\ , but we have the coefficient of the second-order 
term in the hole expansion of lnZg r (Cf fllq)), which is 

U 2 = Z\-\{Z\f (41) 



For N — > oo we may write 



Note that this expression is of order TV, and not of order N 2 as are both terms on its right-hand side. Using ( f40| ) and 
the corresponding expression for Z\ we may equate 

U 2 = y^ e~ /3 ' (ei+ ^ + ^^ ) - - V* e -0*( e <+ e *) (42) 

(id) *d 

U * = -^ E ^^^ - \ E ^^ («) 

and we see that [/ 2 is indeed of order N by virtue of ( J35[) . Note that the terms i = j in the second term of (E3) are not 
compensated by the first term. The second term in ( |43| ) gives the ideal-gas term of the hole system on the two-hole 
level. 

Since we have moments of U2 by our high-temperature expansions, and also the last term in ( p3[ ) is known from 
our one-hole density of states, it is convenient to split 17% into an interacting and an ideal part 

U 2 = C4 nt + Ul d (44) 

with 



dep{e)e- 2pt£ (45) 



2 

N 

~ ~Y 

(Cf. (p4|)). (j>ij must then be determined from U™ 1 . In a continuous version the equation for <j)(e,e') becomes 

jjmt = _Pi N j d£ p( £ ) j de'p(e')e-^ e+E 'U{e,e') (46) 

This relation is not strong enough to yield a unique 4>(e,e'), in the same way as the second virial coefficient of a 
classical gas is not sufficient to determine the interaction potential. The freedom in choice will be reflected upon 
the efficiency of the program to determine the higher-order interactions. We have chosen to have the dependence of 
(j)(e, e') only on the sum variable s + e', and we approximate it by a polynomial: 

_ ( 4? ) 

1 

Equating moments in (f4lf) and 



we have 





0( £ , e ') = E^( £+e ') / 
1 


d in 


k 


u'T) k - 


= Y,^TyJ d ^)J d£ 'P^ + e') k ~ 1+l <i>i 



(48) 



(49) 



Because we are working on a bipartite lattice, all odd moments of L^ 11 * are identically zero. Hence k is even, and 
as also p(e) has only even moments, the combination k — 1 + I must be even and therefore the sums in Eqs. ( |47| ) 
and ( [49|) contain only odd I, The set of equations (^9|) for a finite number of the moments (t/™*), determines an 
equal number of coefficients 4>i- We have computed (U^) , for the square lattice at U = 00, up to k = 12. This 
involves 6 terms (fc = 2, 4, ... , 12) and so we can determine 6 values </>i, ^3, . . . , <pn. In the equations we thus need 
k — 1 + 1 = 12 — 1 + 11 = 22 as the highest moment of p(e), which is just the number of moments we have determined. 
The values of -^ (U^) , are given in Table 0. For the ferromagnetic system (m = 1) these coefficients are zero, as 
the holes do not interact in that case. 

Finally we give the continuum form of the expressions pq) and ( p7[ ) for the grand potential: 

In Zfjp = -pp h N + N j dep{e) ln(l + e~ 0s ) + y N f de f de'p{e)n(e)p{e')n{e')<j)(e, e') (50) 

with 

e = te - [i h + t / de'p(s')4>(s, e')n(e') (51) 
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V. INVERSE SUSCEPTIBILITY 



We return to the uniform susceptibility 



Xfn 



dM 
~~dh 



(52) 



h=0 



with M the total magnetization of the system. As before, we try to find indications of divergences of X fm , which 
should be related to second-order phase transitions between a paramagnetic and a ferromagnetic state. It is usually 
more convenient to express this by stating that the inverse susceptibility must be zero 



Xfa 







and to study 



X f „ 



dh 

dM 



M=0 



or 



PN X 



d(3h 



rm 



sdr 



(53) 



(54) 



(55) 



ro=0 



where m is the magnetization per spin as defined in Sec. [II. In order to find an expression for h, to be able to 
calculate (pq), we construct a generalized (Landau like) free energy 



<p(n B ,fi,m,h) 



1 



\nZv 



(m s 



hmn* 



(56) 



where lnZ gr is given by (22). ip has to be minimized with respect to /i and m at fixed particle density n s and field h, 
to obtain the free energy. Note that this h is not the same field as we used before in Sec. BO. There we interpreted 
h as a field that is felt only by the spins in the background, whereas now we obtain the physical external field that 
would be necessary to yield the given magnetization. Of course, in the case of a finite number of holes (the limit of 
half filling), these fields are the same, as we will see in the resulting expressions. Note also that, due to the definition 
of m as the magnetization per spin, its conjugated variable is hn s , not h. 
We can rewrite (pq ) using (|l^) and ([19]): 



flip = f3e h( n s + /3^h«h 



— In Zt 

N gI 



(57) 



where we can interpret the first term as the contribution of the background of spins, and the other terms as the 
contribution of the holes. Minimization leads to the following equations: 



n h 



dep(e) 



1 



(3h = j3hhi - 



I _|_ e /3(te-/ih) 

cfeM^ln' 
n*dm 



1 



-^(te-Mh) 



with 



hhi = 



de hf 
dm 



The expression for the inverse uniform susceptibility ( p5[ ) then becomes 

d 2 p(e 



i d(3h M 



- / de 



}dm 2 



In 1 



-/3(te-Mh) 



This can be rewritten in terms of p(e), using the legendre transform dlfl) (thus p(e, m) — p(e, /3/ihf)) : 



Wx- 



d(3h 



hf 



rm 



n s dm 



df3h 



hf 



m=0 



n s dm 



m=0 



ds 



d 2 p{e, /3h hi ) 



d{f3h M f 



In 1 



h M =0 



-/3(t£-p h ) 



(58) 
(59) 

(60) 
(61) 

(62) 



Note that m — is equivalent to /ihf = 0, and that, for reasons of symmetry, the first derivative of p with respect to 
hhf vanishes at h^f = 0. 

According to ( p3| ) we want to find values of rih and j3t for which the right-hand side of (p2) is zero, with n^ fixed 
by Eq. (p8[). One can easily verify that, for infinite U, we have Ph^m) = arctan(m), so putting (|62| ) to zero gives 



d 2 p(e,/3h h{ ) 
de 



In (l + e- 13 ^-^) = 1 -n h (63) 

h hi =0 



d(ph h t) 
This equation can be solved by an iterative procedure to calculate the value of /Zh for a given value of [3t. The density 



of states p{e), necessary to calculate rih according to (pq), is determined from its moments as described in Sec. Ill 
and its second derivative is calculated in a similar way. 



To include the interaction described in Sec. IV, one should use the grand potential as given in ( p0[ ) rather than 
the non- interacting hole approximation of (J22|). The final equation, equivalent to (|62j), then involves one extra term 
which contains the second derivative with respect to (3hhi of the interaction (f>. We give a derivation of this equation 
in Appendix B. 

In Fig. we show Curie temperatures for the square-lattice Hubbard model at infinite U, in three different ap- 
proximations: (a) The non-interacting hole approximation, with p determined by interpolation from 8 of its moments 
(of which 4 moments are non-zero); (b) The same but with p determined from 22 (11 non-zero) moments; and (c) 
The interacting-hole approximation, with p determined from 22 moments and <j> from 12 (5 non-zero) interaction 
coefficients. 

One can see that the difference between the 8th- and the 22nd-order non-interacting approximations is small. In 
both approximations, ferromagnetism is stable against paramagnetism for nh < 0.27, at low T. The interaction does 
not change this picture very much. It slightly enhances the stability of the ferromagnetic state, up to rih < 0.29. 
The difference between the non-interacting and the interacting approximations becomes larger with increasing hole 
density, as expected. Numerically, the results agree very well for n^ < 0.06. 

In the next section we will treat the case of finite U. We have been able to calculate 8 moments of the density of 
states in that case, thus we can do an eighth-order approximation at the most. One can then calculate merely two 
coefficients 4>l of the interaction, resulting in an approximation of the interaction which is rather crude. We have seen 
that the picture in the non-interacting hole approximation is qualitatively the same as the one in the interacting-hole 
approximation, in eighth order already. For small rih it agrees rather well also numerically. Therefore, we will not 
include the interaction in the following calculations. 

VI. THE NON-INTERACTING HOLE APPROXIMATION FOR FINITE U 

As we pointed out before, at finite U, excitations in the spin background become possible due to the creation of 
pairs of electrons with opposite spin at the same site. This means that extra empty sites are created, and thus the 
number of empty sites in the system is no longer fixed. Taking U large, however, we can consider the contributions 
to the partition function due to these excitations to be small corrections of the infinite-[7 system, and we can neglect 
the terms that would arise from permanently present electron pairs. To do this, we consider the grand potential of 
the Hubbard model on a square lattice up to the second-order term (taken from Ref. g; note that h is the parameter 
in the Hamiltonian here, not the physical magnetic field we discussed in the previous section): 

_ ^ = i n (i + 2e ^ cosh(/3n) + e 20 ^ u ) + 

4e^(l + e 2 ^-P u ) cosh(/3/i) + 7#7e 2 ^(l - e~W\ 

+ m 2 — — p — - 2 + ■ ■ • ( 64 ) 

In this expression, we will neglect the terms that contain the exponential of —/3U, but we keep terms that are 
proportional to a power of 1/U. This precisely distinguishes the terms that are due to permanent electron pairs, which 
cause an energy @U, from those due to temporary excitations in a system where otherwise no double occupancies 
are present. It is necessary to make this approximation, as the exponential terms can not be treated in this method. 
However, it can be seen easily that these terms are always exponentially smaller than other terms in the expansion, 
and thus that this approximation is justified. 

First we consider the case of half filling, where we have u = U/2: 

3n hf / B „, \ 8e / " 7 / 2 cosh(/3/i) + 7 # 7 (e' 3C/ -l) 
" hf \n 2 + 2e^ u ^cosh((3h))+((3t) 2 [P - puK 5 - + ... (65) 



N 



{2 + 2ef su / 2 cosh{f3h)y 



Here we can neglect all but the most important terms at large U, i.e. we only take the terms containing the highest 
power of e^ u , to get 



pn 



hf 



(3U 



N 



ln{2cosh((3h)) + {f3t) 



By definition, this expression must be equal to 



pv 



{flU) (cosh(/3/i)) 
^ ln^N, so using the definition ([12]) for ehf we get 



(66) 



\ 



/?e hf = In (2 cosh((3h)) + {(3tf 



((3U) (cosh(/3h)y 



(67) 



and we see that this is indeed a correction of order jj on Eq. (|T|). Note that we obtain the same result if we first 
omit the e~° u terms in (64), and only then substitute U/2 for /i. This once more supports our statement that these 
terms may be neglected. 

Off half filling, we have to rewrite ( p4| ) (without the e~^ u terms) in terms of the effective chemical potential /^h for 
the holes, as defined by Eq. (|l6|), but now containing the corrected ehf as given by (J67|)- For simplicity, we do this in a 
few steps. First, we substitute the chemical potential for the holes without the correction terms, as in Sec. [II. Then 
we expand the logarithm and the numerators with respect to the exponential of this chemical potential. Finally, we 
include the corrected /ih by expanding the exponentials with respect to the correction terms. Thus, we obtain for the 
grand potential 



£} = -/3 M h + e^ 



1 + {ptf 



2- 



(0U) (cosh(ph)Y 



(68) 



The coefficient of e^ Mh in this expression again determines the moments of the distribution p(e, /3h), as described in 



Sec. III. Of course these are now functions of f3U. In Table III we give the moments that we have been able to derive 
from the series expansion data, for the square lattice, at h = 0. Note that the moments for h = oo are the same as in 
the case of infinite U (Table [fl), because U has no significance in a system where all spins point in the same direction. 
We can now apply the method described in Sec. |V| to calculate Curie temperatures for finite U. One has to realize, 
though, that at half filling the inverse susceptibility depends on the temperature, which was not the case for infinite U. 
Due to the excitations we get corrections of the type (3t 2 /U , thus we still have a series expansion in the parameter fit. 
The coefficients in this expansion are suppressed by large factors ftU, however, and the range of convergence of the 
expansion is (3t < 30 or further, depending on the value of j3U . Thus, we may hope that convergence is good enough 
in the region where we expect to find solutions of (p3). We give the full expression for the inverse susceptibility at 
half filling, for the square lattice and up to the (/3t)°terms: 



WXh? = 



d(3h 



hf 



dm 



= 1 



4(/3t) 2 8 (-2 + 0U) (f3t) 4 



(/?[/) 0917)3 

(1131 - 648/3C/ + 32(/3£f) 2 ) {(3tf 

3(/?£/) 5 + 

(-9129 + 6296/JC/ - 1132(/3C/) 2 + A{f3Uf) {fitf 

wy 



(69) 



We have checked that (|69j) does not become zero for any value of (3t and (3U. Therefore we expect no transition from 
a paramagnetic to a ferromagnetic state in the half-filled system. Thus, we only have to consider the second factor 
on the right-hand side of (B2), which vanishes at 



de 



d 2 p(s,(3h) 



d{phy 



In 1 



-0(te-fi h ) 



h=0 



Xhf 






(70) 



We show the results for the square and the simple cubic lattices in the next section. 



VII. MAGNETIC-PHASE DIAGRAM 

We have used the theory described above to calculate Curie temperatures for the square and simple cubic lattices. 
For both lattices, we find a surface of Curie temperatures in the n^-jj-T diagram. In Figs, ra and o we display these 
results in various ways. 
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In Figs. |p, and ||a, contours of fixed Curie temperature are plotted in the n^-fj plane. In the range of temperatures 
up to about ^j- = 0.20 we find a curve enclosing a region of ferromagnetism. For -p > 0.07 these curves are closed 
and lie away from the 4? = axis. Thus, at given density n^ and temperature Tq, one has to go to finite U to find 
a transition. In other words: allowing for excitations in the spin background enhances the ferromagnetic behavior. 
Furthermore, curves are generally not enclosed by all contours at lower temperatures. This would imply that, at given 
nh and 4j, one would find a para-ferro transition when lowering the temperature from a region of high temperature, 
but also when letting it increase from zero. This reentering of a paramagnetic phase at low temperatures does not 
seem to be physical. It is in fact an artefact of this method, due to convergence problems at very low temperatures. 
One can understand this by looking at the expression (pjj). If the highest-order term becomes of order one, the series 
is clearly too short and does not converge properly. This means that the results become unreliable for jj > ^ in 
the case of the square lattice, and ^>^on the simple cubic lattice. For a few curves we have indicated this by a 
dashed line. As the approximations are better for higher temperatures, we assume that the actual curve at Tq = 
(for which we can only perform a calculation at infinite U) should enclose all curves shown. 

In Figs, pb and fib, we show Curie temperatures in contours of fixed ij. Again we see the non-physical behavior 
of curves being closed at the low- temperature side, for almost all values of jj. Fig. [?] shows Curie temperatures at 
fixed rih = 0.09, for the simple cubic lattice, indicated by the dotted-dashed lines in Figs. Ha and Kp. The dotted line 
in Fig. [] indicates the region where the series expansion becomes unreliable, according to the arguments presented 
above. 

There is one other point we want to mention here. As we have stated in the introduction, we have also constructed 
the staggered susceptibility by replacing the magnetic field ft by a staggered field h s . Although it is much more 
complicated to calculate the high-temperature expansions for that case, as the number of terms involved increases 
significantly, it is not difficult to obtain expressions for the staggered susceptibility, both at half filling and in the 
one-hole approximation, for h s = 0. Thus, one may think that it is possible to obtain similar results for the transition 
between a paramagnetic and an antiferromagnetic state, and conclude which transition occurs first. When putting 
the inverse staggered susceptibility at half filling (the equivalent of (p9) for the antiferromagnetic system) to zero, 
one finds solutions for all values of the parameter (3U. This means that the staggered susceptibility of the half-filled 
system diverges at a finite temperature. Apparently, the para-antiferro transition is driven by the background itself, 
and may be disturbed by a finite hole density. In our formulation, however, it is the holes that drive the system into 
an ordered state, and the background only indirectly contributes to the transition via its interaction with the holes. 
This formulation is clearly not suitable to describe the transition to an antiferromagnetic state. Therefore we only 
briefly indicate what we expect for the para-antiferro transition. 

In Fig. H we plot Neel temperatures for the simple cubic lattice at half filling, in approximations to different orders 
in the parameter f3t. We see that the convergence of the series expansion is very good for large U. A transition from 
a paramagnetic to an antiferromagnetic phase is expected for all values of U. It is at Tn = for infinite U, and 
at increasing temperatures with decreasing U. For finite hole densities we expect the transition to occur at lower 
temperatures, and at some point cross the para-ferro transition. 

VIII. DISCUSSION AND CONCLUSIONS 

We have calculated Curie temperatures for the large-?/ Hubbard model on the square and simple cubic lattices, 
by means of a new extrapolation method to extract information on low-temperature behavior from high-temperature 
series expansions. We find a region of ferromagnetic behavior in the magnetic-phase diagram, near half filling. 

Comparing previous results for the simple cubic lattice, as depicted in Fig. El to our current results, shown in Fig. 0, 
we see that we now find a Curie temperature that is an order of magnitude smaller than before. Furthermore, as we 
have checked in the case of infinite U, subsequent approximations in the current method do give consistent results, 
instead of alternatingly producing Curie temperatures or not. These convergence problems in the primitive series 
expansions are likely due to the fermi degeneracy of the electron gas. At j3t ~ 1, the wavelength of the electrons 
becomes equal to the lattice distance, causing this degeneracy and divergences to be present. When applying a 
straightforward extrapolation technique, one can not account for this degeneracy, leading to results that are erroneous 
for j3t > 1. In our approximation, using a density of states for holes, we take the fermi degeneracy into account, and 
therefore we are able to proceed to lower temperatures. We are confident that our present results do not suffer from 
the above-mentioned convergence problems. 

As we show in Fig. 0, the difference between approximations to different orders in the parameter f3t is rather 
small, and adding the interaction also does not change the result considerably. Thus we believe the eighth-order non- 
interacting hole approximation to be sufficient to describe the qualitative behavior, and to obtain a good indication 
for numerical values. We may add that, as a check, we have compared the free energy from calculations by this 

If 



method to results following directly from the series expansions, at /3t < 0.5, where the expansions are almost exact, 
and that these results agree very well. 

Our method works only for large U, low hole density (rih < 0.2), and, depending on the value of U, sufficiently 
high temperature. This is clear from Figs. f^-fl where we see that the results are unreliable for ir > ^. We believe, 
however, that our method gives a correct description for the tendencies in the half-filled system at infinite U, and for 
the qualitative behavior up to n^ rs 0.2. 

There are, however, some important limitations to this method, due to which we are not able to predict a ferro- 
magnetic state with certainty. ■— ■ 

As we know from a theorem by Ghosnl3, similar to the Mermin- Wagner theorem for the Heisenberg model, the 
Hubbard model does not have long-range ordering in two dimensions for finite temperatures. Thus, we must expect 
a ferromagnetic phase in the two-dimensional case to be of the Kosterlitz-Thouless type. Our method is essentially 
based on short-range information from the high-temperature expansion (which is obtained via calculations on small 
systems). It gives similar results for the square and the simple cubic lattices, as can be seen in Figs. H and o, and we 
can not distinguish between different kinds of phases occurring. 

Also, the method currently fails to describe the case of a para-antiferro transition, due to the fact that a divergent 
background is not treated correctly. We can therefore calculate only possible second-order phase transitions between 
a paramagnetic and a ferromagnetic phase, for the case of finite hole density. At half filling, we do find a finite 
Neel temperature for any finite value of the parameter (3U (see Fig. ||). This implies that, near half filling, there is 
a transition from a paramagnetic to an antiferromagnetic state at a higher temperature than the calculated para- 
ferro transition. Thus, the para-ferro transition can not occur, and one must study the antiferro-ferro transition to 
determine the ground-state behavior. 

Finally, due to the thermodynamic approach in which all possible states are taken into account, our method can 
not distinguish special states that may start to dominate the system at low temperatures. Such states, if any, are not 
recognized by the high-temperature expansion. An example of this is the fact that it fails to reflect the influence of 
m = 1 states in an m = system. 

We can compare our results to the work of Putikka et al.EJ, who calculate series expansions similar to those used 
by us, for the related t—J model, and extrapolate to low temperatures by means of Pade approximants. For J > 0, 
in the limit of small J, the t—J model is equivalent to the large-[/ Hubbard model. They find a region of weak 
ferromagnetism (i.e. the spins are not fully aligned) for small positive J, at hole density rih < 0.28 ±0.05, which is in 
good agreement with our results. 

It is also encouraging to note that some of our results are in reasonable quantitative agreement with results using an 
approximation-of an entirely different nature. By means of the slave-boson mean-field approach (at T = 0), Dcntcnccr 
and Blaauboerli3 find a critical hole density njj = 1/3 for ferromagnetism to occur at U = oo, in agreement with the 
values 0.27-0.29 found here (see Fig. 0). They also find that the value of U/t above which ferromagnetism can occur, 
is U/t = 20 (at rih = 0.17), whereas one may extrapolate the results of our Fig. ||a to T = to find U/t = 15 (at 
n h = 0.15). 

Von dcr Linden and EdwardsEj use a variational approach to find a ferromagnetic region in the T — phase 
diagram of the square-lattice Hubbard model. They rigorously conclude that the state of complete spin alignment is 
unstable when n^ > 0.29, for all U, and when U/t < 42, for all nt- The latter value is significantly higher than the 
value above which we find ferromagnetism, but we assume that that is due to the fact that they consider only strong 
ferromagnetism (full alignment of the.spins), whereas our method may also include weak ferromagnetism. 

Also the results of Barbieri et al.E-3, who consider systems with a large (but finite) number of holes, support the 
existence of ferromagnetic behavior. 

A final comparison that can be made is for the relation between the Neel temperature and U/t in the half-filled 
system. From Fig. one can calculate that the para-antiferro transition occurs for kT^ m 3.85t/U. The large-I/ 
Hubbard model at half filling is known to be equivalent to an antiferromagnetic Heisenberg spin model, for which 
estimates of the values of the critical temperature are given in Ref. H3. According to the results mentioned there, the 
relation would be kT^ ss 3.80i/C/, which is in very good agreement. 

This research was supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially 
supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). 
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APPENDIX A: ENUMERATION OF PATHS 

In this appendix we describe an efficient way to calculate the moments of the density of states, for the case of 
infinite U, by which we have calculated 22 of these moments for the square lattice, and 16 for the simple cubic lattice. 
We start from Eq. (£0|), which we expand in terms of the parameter (3t: 

-^ = Y,W -^-M n (m) (A.l) 

71 = 

with the moments of the density of states defined as 

M n {m)= J dep{e,m)e n (A.2) 



We can write the partition function for one hole according to its definition (Cf. (|10[)) also as 

Z\ = e [N - 1)l3 ^Z N ^ (A.3) 

N^ra)) £ <*,a,(m)|e-^|t,a,(m)> (A.4) 

\i,ai(m)) 

where the summation is over all states \i,ai(m)) with a hole at site i and with a spin background ai(m) such that 
the magnetization per spin is indeed to. N^ denotes the number of electrons with spin up, which depends on to, and 

the factor I „ , , ) is the total number of possible background configurations given the location of the hole, which 

accounts for the spin degrees of freedom. In the thermodynamic limit, this factor is exactly equal to the exponential 



summation over i gives a trivial (translational) factor N, and we can expand the exponential in powers of f3t to obtain 



factor in (A.3), as one easily checks by applying Stirling's formula for the binomial, and with ( |l3| ) for eht- The 

(translation 

Zl(N-l' 

V IV '' |a(m)) n 

where 

A n (a(m)) = {a(m)\C^] |a(m)> (A.6) 

is the number of w alks of len gth n in the configuration space, that restore the spin background a(m) to its original 



state. Comparing (A.l) and (A.J;) we see that 



M »( m ) = (tf t J)) £ M<x(m)) (A.7) 

\a(m)) 

Thus M n (m) is precisely the sum over all possible closed walks w n of length n, summing the fraction of spin back- 
grounds that is restored by w n . Such a walk induces a permutation V(w n ) of the background spins, which can be 
written as a product of disjunct cyclic permutations Vi{w n ) with length \Vi{w n )\ > 1. In order to restore the spin 
background a(m), the direction of the spin on each site must remain unchanged, when applying Vi(w n ). Thus, all 
spins that are interchanged by this permutation must point in the same direction. As the number of spins involved 
is negligible compared to the total number of spins, we may approximate that the probability to find an individual 
spin pointing up or down is given by -4^ anc i 1^S ) respectively. Hence the fraction of backgrounds in which the 

alignment of the spins remains unchanged under the permutation Vi{w n ) is (— j 2 ) + (^r 1 ) j where / = \Vi(w n )\ is 
the number of spins involved in the permutation. Thus, we can calculate M n as 

M„ (m ,^n((_t-) + (v) ) (A.S, 

For the actual evaluation of this expression we proved an elegant theorem that enables us to significantly extend 
earlier calculations of the moments to n = 22. Defining a retracing sequence as two subsequent steps of the hole in 
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opposite directions (thus after two steps the hole is back in its previous position; note that the last and first steps of 
a closed walk are considered to be subsequent as well), one can make a distinction between reducible and irreducible 
closed walks: an irreducible walk does not contain any retracing sequence, whereas a reducible walk does. A reducible 
walk can be made irreducible by repeatedly removing its retracing sequences; the result is called the irreducible part 
of the walk. Note that a retracing sequence does not permute spins, so the irreducible part of a walk induces the 
same permutation of the spins as the walk itself. Thus, it is sufficient to study only irreducible walks if one knows of 
how many reducible walks of a given length it is the irreducible part. We proved the following formula: the number 
of closed walks of length I + In on a hypercubical lattice with coordination number z, that have a given irreducible 
part of length I > 0, is 

N il (l,n) = (z-ir( l+ n 2n ) (A.9) 



This greatly facilitates the calculation of ( |A.8| ). 

APPENDIX B: INVERSE SUSCEPTIBILITY IN THE INTERACTING-HOLE APPROXIMATION 

In this appendix we g ive the formula for the inverse susceptibility in the interacting-holc approximation, using the 
theory given in Sec. [fyL We start from Eq. (|5q), which has to be differentiated with respect to m in order to get the 
equivalent of (|j9|), with ( |50| ) for lnZ gr : 

f3h = (3h M + n h ^- fde^-ln(l + e-^)+ f dep(e)n(i)- 90§ 



n s dm J n s dm J n s dm 

<3t fde f de'^P-n(i)p(e')n(e')<f>(e,e') 
J J n s dm 

/3t [ de J de'p(s)^P-p(s')n(e'U(e,s') 
J J n s dm 

| fde fde'p(e)n(s)p(e')n(e')^^- (B.l) 



where n^ is given by 



n, = - J dep(e)n(i)-^- + fit J de J de'p{e)n[e)p{e')^l<l>{e,e') (B.2) 



This may look awkward, but if we look at the derivatives of e (see Eq. (|51|)) we see that many of these terms cancel. 



Let us first look at the expression (B.2) for the hole density. As we are working at fixed hole density, derivatives of 
the fermi factor do not play a role in these equations, and they vanish. We need the derivative of e with respect to 

dl3£ =-l + (3t J de'p(e')(p{e,e')^^- (B.3) 



dfinh J " " ' " d/3p h 

and so we see that indeed there is a cancellation of terms, leaving us with the relation 



n h = / dep(e)n(e) (B.4) 



Then, we rewrite the expression for the magnetic field with 



«' ' 9 ^ h +f* [de'^fin(e'me,e') + (3t [ de>p(e')^^e>) 



dm n s dm J n s dm ' J n s dm 

+ 0t fde'pie'Me') 9 ^^ (B.5) 

./ n s dm 



Using this expression it is straightforward to check that (B.l) reduces to 



(3h = (3h h{ ~ I de^^-ln(l + e- 0S ) + $ [ de [ de'p(e)n(e)p{e')n(e')^^- (B.6) 

n s om 2 J J n s om 
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In order to derive the inverse susceptibility from this expression, we have to take the derivative with respect to n s m 
again, and put m = 0. For reasons of symmetry it is easy to show that the firs t de rivatives with respect to m of 
all functions appearing in the integrals vanish at m — 0. Thus, in the terms in (B.6) we only have to consider the 
derivatives of the functions that have been differentiated once already: 



PNxi' = 



d(5h 



hf 



n s dm 



m=0 



ds 



ftt_ 
~2 



d 2 p(e) 



n 2 dm 2 



ln(l 



,-Pe\ 



ro=0 



ds / ds' p(s)n(s)p(s')n(s') 



d 2 <P(s,s') 



nidm 2 



(B.7) 



ro=0 



This can again be expressed in terms of p(s) (note that also </> is being legendre transformed): 



PNxil = 



d[3h 



hf 



n s dm 



1- 



d(ih 



iif 



n s dm 



f ' ■'■■ 



^m$^ 



,-0S\ 



d 2 

ds' p(s)n(s)p(s')n(s') 



£,£' 



d(f3h) 2 



h=0 



(B- 



which is the modification of (162) for interacting holes 



1 Y. Nagaoka, Phys. Rev. 147, 392 (1966). 

2 For a summary, see e.g. E. Muller-Hartmann, Th. Hanish, and R. Hirsch, Physica B 186—188, 834 (1993), and references 
therein. 

3 D.F.B. ten Haaf and J.M.J, van Leeuwen, Phys. Rev. B 46, 6313 (1992). 

4 J.A. Henderson, J. Oitmaa, and M.C.B. Ashley, Phys. Rev. B 46, 6328 (1992). 

5 K. Kubo, Prog. Theor. Phys. 64, 758 (1980); K. Kubo and M. Tada, ibid. 69, 1345 (1983) and 71, 479 (1984); K.K. Pan 
and Y.L. Wang, Phys. Rev. B 43, 3706 (1991) and J. Appl. Phys. 69, 4656 (1991). 

6 W.F. Brinkman and T.M. Rice, Phys. Rev. B 2, 1324 (1970). 

7 S.V. Meshkov and D.V. Berkov, Phys. Lett. A 170, 448 (1992). 

8 See e.g. A Modern Course in Statistical Physics by L.E. Reichl, Edward Arnold (Publishers) Ltd 1980. 

9 See e.g. A Method for Studying Model Hamiltonians by N.N. Bogolubov, Jr., Pergamon Press 1972. 
D.K. Ghosh, Phys. Rev. Lett. 27, 1584 (1971). 

W.O. Putikka, M.U. Luchini, and M. Ogata, Phys. Rev. Lett. 69, 2288 ( 1992). 

P.J.H. Denteneer and M. Blaauboer, submitted for publication (preprint cond-mat 9406093). 

W. von der Linden and D.M. Edwards, J. Phys.: Condens. Matter 3, 4917 (1991). 

A. Barbieri, J.A. Riera, and A. P. Young, Phys. Rev. B 41, 11697 (1990). 

G.S. Rushbrooke, G.A. Baker, Jr., and P.J. Wood, in Phase Transitions and Critical Phenomena, Volume 3, edited by C. 

Domb and M.S. Green; pp. 305-307. 



15 



TABLES 



n 


m = 


m — 1 





1 


1 


2 


4 


4 


4 


30 


36 


6 


269.5 


400 


8 


2641.75 


4900 


10 


27279.94 


63504 


12 


291718.97 


853776 


14 


3199250.29 


11778624 


16 


35766660.22 


165636900 


18 


405989247.14 


2363904400 


20 


4665921461.01 


34134779536 


22 


54182396281.84 


497634306624 



TABLE pi The first 22 moments of p(e) for to — and to = 1, for the square lattice (odd moments vanish). 



k 


m = 


m — 1 











2 








4 


-0.333333333333 





6 


-0.5875 





8 


-0.424851190476 





10 


-0.170455970293 





12 


-0.041839059880 






TABLE [f| Values of -^ (t/™*) fc for to = and to = 1, for the square lattice with U = oo. 



f dep(e)e 



1 
2 f 2 



/scrj 



24 | 2 £_ 

^U /3C7 



(0uy 



+ 



720 
40320 



1440 
' 10567 



59 

48 -,T 



WW 

93 



swuy 2 



271 



6(/3E7) 
1459 



T + 



576/31/ 320(/3!7) 2 



2(/3!7)' 



+ 



127 
WW? 



377 , 4531 , 4043 
32(/3C/) a "^ 960C/) 4 "^ 8(/3C/) & 



28837 , 78593 
8(J3U)° ~*~ 8(0U)'' 



TABLE III. Moments of the density of states for the square lattice (odd moments vanish), at h = 0. 
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FIGURE CAPTIONS 

FIG. [0. The inverse ferromagnetic susceptibility as a function of the parameter j3t, for the Hubbard model on a 
simple cubic lattice, with infinite U and particle density n — 0.9. Approximations up to order 2, 4, 6, 8 in (3t obtained 
by means of the cluster expansion method. 

FIG. 0. Neel and Curie temperatures, as a function of the particle density, for the Hubbard model on a simple 
cubic lattice, at constant t/U . 

FIG. 0a. The density of states for a paramagnetic system on a square lattice, using up to the number of moments 
indicated. 



FIG. 0b. The density of states in the ferromagnetic regime on a square lattice. Exact result, and approximations 
using up to the number of moments indicated. 

FIG. [|. Curie temperatures (contours of zero inverse ferromagnetic susceptibility) for the square lattice at infinite U. 

(a) Non-interacting hole approximation, 8th order; 

(b) Non-interacting hole approximation, 22nd order; 

(c) Interacting-hole approximation, 22nd order. 

FIG. ||. Magnetic-phase diagram for the square lattice. 

(a) Contours of fixed Curie temperature, with kTc/t — 0.03, 0.04, . . . , 0.19 (increment 0.01). 

(b) Curie temperature at fixed t/U = 0, 0.005, . . . , 0.055 (increment 0.005). 

FIG. g. Phase diagram for the simple cubic lattice. 

(a) Contours of fixed Curie temperature, with kTc/t = 0.03, 0.04, . . . , 0.14 (increment 0.01). 

(b) Curie temperature at fixed t/U = 0, 0.002, . . . , 0.022 (increment 0.002). 

FIG. pi. Curie temperature for the simple cubic lattice, at n^ — 0.09. The dashed part of the curve is unreliable, 
due to lack of convergence (indicated by the dotted line). 

FIG. g Neel temperature for the simple cubic lattice at half filling. Approximations to different orders in (3t, as 
indicated. 
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